#Load data
rm(list=ls())
library(readxl)
mydata<-read_excel("Data_summary_NP.xlsx",sheet = "Sheet2")
head(mydata)
## # A tibble: 6 × 33
## Forest N P Treatment Total_C Total_N TC_TN WSOC NH4 NO3 availP
## <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 CK 33.3 2.24 14.9 55.6 3.53 13.4 0.309
## 2 1 1 1 CK 29.2 1.72 17.0 49.9 2.80 10.2 0.275
## 3 1 1 1 CK 34.8 2.22 15.7 62.6 1.78 16.6 0.710
## 4 1 1 1 CK 37.2 2.21 16.8 107. 1.97 15 0.615
## 5 1 1 1 CK 29.3 1.83 16.0 46.0 2.99 14.2 0.553
## 6 1 2 1 N 22.9 1.61 14.2 39.6 1.39 11.0 1.11
## # ℹ 22 more variables: availS <dbl>, pH <dbl>, EC <dbl>, MBC <dbl>, MBN <dbl>,
## # MBC_MBN <dbl>, ARS <dbl>, SARSmbc <dbl>, SARSt <dbl>, BG <dbl>, SBG <dbl>,
## # NAG <dbl>, SNAG <dbl>, XYL <dbl>, SXYL <dbl>, LAP <dbl>, SLAP <dbl>,
## # ACP <dbl>, SACP <dbl>, logS_C <dbl>, logS_N <dbl>, logS_P <dbl>
str(mydata)
## tibble [60 × 33] (S3: tbl_df/tbl/data.frame)
## $ Forest : num [1:60] 1 1 1 1 1 1 1 1 1 1 ...
## $ N : num [1:60] 1 1 1 1 1 2 2 2 2 2 ...
## $ P : num [1:60] 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment: chr [1:60] "CK" "CK" "CK" "CK" ...
## $ Total_C : num [1:60] 33.3 29.2 34.8 37.2 29.3 ...
## $ Total_N : num [1:60] 2.24 1.72 2.22 2.21 1.83 ...
## $ TC_TN : num [1:60] 14.9 17 15.7 16.8 16 ...
## $ WSOC : num [1:60] 55.6 49.9 62.6 107.3 46 ...
## $ NH4 : num [1:60] 3.53 2.8 1.78 1.97 2.99 ...
## $ NO3 : num [1:60] 13.3 10.2 16.6 15 14.2 ...
## $ availP : num [1:60] 0.309 0.275 0.71 0.615 0.553 ...
## $ availS : num [1:60] 48.6 41.3 32.9 29.6 33.3 ...
## $ pH : num [1:60] 3.7 3.66 3.56 3.58 3.64 3.66 3.48 3.55 3.68 3.5 ...
## $ EC : num [1:60] 119 129 151 146 154 ...
## $ MBC : num [1:60] 363 356 399 349 343 ...
## $ MBN : num [1:60] 70.4 72.7 80.9 50.9 65.8 ...
## $ MBC_MBN : num [1:60] 6.02 5.71 5.76 8.01 6.09 ...
## $ ARS : num [1:60] 49.1 46.6 52.5 53.8 52.4 ...
## $ SARSmbc : num [1:60] 0.135 0.131 0.132 0.154 0.153 ...
## $ SARSt : num [1:60] 1.47 1.59 1.51 1.45 1.79 ...
## $ BG : num [1:60] 244 278 312 255 326 ...
## $ SBG : num [1:60] 0.67 0.781 0.781 0.731 0.949 ...
## $ NAG : num [1:60] 216 178 372 172 239 ...
## $ SNAG : num [1:60] 0.595 0.501 0.932 0.492 0.696 ...
## $ XYL : num [1:60] 300 310 306 337 323 ...
## $ SXYL : num [1:60] 0.826 0.872 0.766 0.965 0.94 ...
## $ LAP : num [1:60] 21.6 24.8 24.6 18.8 18 ...
## $ SLAP : num [1:60] 0.0593 0.0695 0.0615 0.0538 0.0524 ...
## $ ACP : num [1:60] 18969 9279 6488 9556 10707 ...
## $ SACP : num [1:60] 52.2 26.1 16.2 27.4 31.2 ...
## $ logS_C : num [1:60] 0.618 0.602 0.616 0.624 0.611 ...
## $ logS_N : num [1:60] 0.712 0.723 0.662 0.759 0.713 ...
## $ logS_P : num [1:60] 0.395 0.42 0.451 0.435 0.427 ...
mydata$Forest<-factor(mydata$Forest,levels = c(1,2,3),labels = c("PF","RF","DF"))
mydata$Forest2<-factor(mydata$Forest,levels = c("PF","RF","DF"),labels = c("Primary Forest","Rehabilitated Forest","Disturbed Forest"))
mydata$N<-factor(mydata$N,levels = c(1,2),labels = c("N0","N150"))
mydata$P<-factor(mydata$P,levels = c(1,2),labels = c("P0","P150"))
mydata$Treatment<-factor(mydata$Treatment,levels = c("CK","N","P","NP"),labels = c("CK","N","P","NP"))
#data transform for regression
mydata.log<-data.frame(mydata[1:4],log(mydata[,c(5:12)]+1),mydata[13],
log(mydata[,c(14:30)]+1),mydata[31:34])
#data transform for statistical analysis
mydata.log2<-data.frame(mydata[1:4],log(mydata[,c(5:12)]),mydata[13],
log(mydata[,c(14:30)]),mydata[31:34])
#plot
library(vegan)
library(ggpubr)
library(readxl)
library(ggplot2)
#Fig. 1a
#PCA
library(vegan)
set.seed(123)
adonis<-adonis2(mydata[,5:17]~Forest*N*P,data = mydata,permutations = 999,method = "bray")
adonis
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = mydata[, 5:17] ~ Forest * N * P, data = mydata, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Forest 2 1.5673 0.45191 101.6490 0.001 ***
## N 1 0.0264 0.00760 3.4199 0.035 *
## P 1 1.7250 0.49739 223.7575 0.001 ***
## Forest:N 2 0.0258 0.00744 1.6739 0.159
## Forest:P 2 -0.2705 -0.07800 -17.5457 1.000
## N:P 1 0.0151 0.00436 1.9633 0.151
## Forest:N:P 2 0.0090 0.00259 0.5834 0.651
## Residual 48 0.3700 0.10670
## Total 59 3.4681 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
data.ord<-mydata
colnames(data.ord)
## [1] "Forest" "N" "P" "Treatment" "Total_C" "Total_N"
## [7] "TC_TN" "WSOC" "NH4" "NO3" "availP" "availS"
## [13] "pH" "EC" "MBC" "MBN" "MBC_MBN" "ARS"
## [19] "SARSmbc" "SARSt" "BG" "SBG" "NAG" "SNAG"
## [25] "XYL" "SXYL" "LAP" "SLAP" "ACP" "SACP"
## [31] "logS_C" "logS_N" "logS_P" "Forest2"
pca<-prcomp(data.ord[,c(5:17)],scale. = TRUE)
pca
## Standard deviations (1, .., p=13):
## [1] 2.61106181 1.57448982 1.06089326 0.92335728 0.80663187 0.56308332
## [7] 0.54809885 0.46327968 0.32990260 0.28215176 0.21973926 0.05681361
## [13] 0.05038036
##
## Rotation (n x k) = (13 x 13):
## PC1 PC2 PC3 PC4 PC5
## Total_C -0.36439584 0.03837000 0.137217479 -0.20596512 0.05867890
## Total_N -0.37210753 -0.02244927 -0.014407924 0.01181656 0.06315191
## TC_TN -0.11432172 0.20222533 0.574874493 -0.70339382 -0.02857393
## WSOC -0.23280945 -0.41126889 -0.112238039 -0.24623510 -0.01600838
## NH4 -0.23525969 -0.20310170 -0.264447430 -0.10875217 -0.73561753
## NO3 -0.31843461 0.23047746 -0.003701715 0.19693792 0.10081928
## availP -0.02944763 -0.53330539 0.127331621 -0.01160790 0.46303372
## availS -0.17094228 0.48920314 -0.018176389 0.07301420 -0.15682820
## pH 0.30569171 -0.25344614 -0.006592140 -0.22109512 -0.25774976
## EC -0.36528912 0.04065440 0.011717626 0.17741966 0.13229570
## MBC -0.35840990 -0.16210917 -0.018562525 0.04648275 -0.03071429
## MBN -0.32891341 -0.23353058 0.167594344 0.18421945 -0.10376330
## MBC_MBN -0.10545405 0.15224465 -0.723037106 -0.47907047 0.32489325
## PC6 PC7 PC8 PC9 PC10
## Total_C 0.03606198 -0.118315994 -0.11137249 -0.35198276 0.341692996
## Total_N 0.08152112 -0.201517143 -0.11313135 -0.38670423 0.456420136
## TC_TN -0.08308426 0.173696760 0.07011910 0.08706354 -0.155272968
## WSOC -0.17190205 -0.556128055 -0.17551868 0.57618894 0.036580507
## NH4 -0.10139456 0.486205945 -0.13915271 0.03073294 0.103179534
## NO3 -0.38311374 -0.003469945 -0.54571942 -0.13189383 -0.530475844
## availP 0.33208727 0.479787145 -0.34354163 0.09513967 0.007172691
## availS 0.70067226 -0.067396593 -0.26346282 0.36589470 0.004477665
## pH 0.30924600 -0.301126835 -0.34837681 -0.43738973 -0.327174711
## EC -0.09232694 0.150071956 -0.01889438 0.09023699 0.010526538
## MBC 0.19225484 -0.083195462 0.37711873 -0.13473744 -0.313920335
## MBN 0.22029948 -0.051106031 0.38604978 -0.03722634 -0.347135430
## MBC_MBN 0.08757186 0.109839186 0.14659785 -0.08204118 -0.176503570
## PC11 PC12 PC13
## Total_C -0.02756986 0.203396841 -0.700407468
## Total_N -0.12602835 -0.272794833 0.590324834
## TC_TN 0.02579549 -0.065523555 0.203314281
## WSOC -0.01001192 0.003265111 -0.009696621
## NH4 -0.05331246 0.005206253 0.003089812
## NO3 -0.22978057 -0.016179737 0.007576697
## availP -0.12130498 0.037661382 0.015332711
## availS -0.04213588 0.024395028 -0.006031454
## pH 0.34532981 -0.027943764 0.029878391
## EC 0.87957675 -0.008370397 0.049114101
## MBC -0.09264389 0.687755329 0.239628869
## MBN -0.09401326 -0.613832381 -0.238311408
## MBC_MBN 0.01140767 -0.164188364 -0.040713501
summary(pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.6111 1.5745 1.06089 0.92336 0.80663 0.56308 0.54810
## Proportion of Variance 0.5244 0.1907 0.08658 0.06558 0.05005 0.02439 0.02311
## Cumulative Proportion 0.5244 0.7151 0.80170 0.86729 0.91734 0.94173 0.96484
## PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 0.46328 0.32990 0.28215 0.21974 0.05681 0.05038
## Proportion of Variance 0.01651 0.00837 0.00612 0.00371 0.00025 0.00020
## Cumulative Proportion 0.98135 0.98972 0.99584 0.99956 0.99980 1.00000
library(ggbiplot)
ggbiplot(pca,obs.scale = 1,var.scale = 1,ellipse = TRUE) +
scale_color_manual(name = "Treatment",values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
scale_shape_discrete(name="Forest")+
geom_point(aes(shape=mydata$Forest,colour=mydata$Treatment),size=5)+theme_bw()+
theme(legend.direction='vertical',legend.position='right')+
theme(axis.text.x = element_text(size = 15,color="black"))+theme(axis.text.y = element_text(size = 15,color="black"))+
theme(axis.title = element_text(size = 15,color="black"))+
annotate("text",x=1.25,y=-2.5,label="Permanova analysis:\nForest: R2 = 0.45;p = 0.001\nN: R2 = 0.0076 p = 0.035\nP: R2 = 0.50; p = 0.001",color="black",size=3)

#Fig. 1b
library(ggpubr)
label = "p.signif"
method = "wilcox.test"
my_comparison<-list(c("CK","N"),c("CK","P"),c("CK","NP"),c("P","NP"))
size=17
g1<-ggboxplot(mydata,x='Treatment',y='availS',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Extractable S (mg kg-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g1

#ggsave("plot2/Fig.1b_availS.pdf",plot = last_plot(),units = "in",height = 5,width = 7,dpi = 600)
#Fig. 2
#ARS#
g3<-ggboxplot(mydata,x='Treatment',y='ARS',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Absolute ARS activity\n(nmol/g/h)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g3

#"BG"
g4<-ggboxplot(mydata,x='Treatment',y='BG',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Absolute BG activity\n(nmol/g/h)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g4

#"XYL"
g5<-ggboxplot(mydata,x='Treatment',y='XYL',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Absolute XYL activity\n(nmol/g/h)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g5

#"NAG"
g6<-ggboxplot(mydata.log,x='Treatment',y='NAG',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Absolute NAG activity\nLog(nmol/g/h)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g6

#"LAP"
g7<-ggboxplot(mydata,x='Treatment',y='LAP',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Absolute LAP activity\n(nmol/g/h)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g7

#"ACP"
g8<-ggboxplot(mydata.log,x='Treatment',y='ACP',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Absolute AP activity\nLog(nmol/g/h)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g8

ggarrange(g3,g4,g5,g6,g7,g8,ncol = 2,nrow = 3,widths=c(1,1),heights=c(1,1),
common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)","(e)","(f)"))

#ggsave("plot2/Fig.3_enzymeActivity2.pdf",plot = last_plot(),units = "in",height = 15,width = 14,dpi = 600)
#Fig. 3
g9<-ggboxplot(mydata.log,x='Treatment',y='SARSmbc',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific ARS activity\nLog(nmol ug MBC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g9

g10<-ggboxplot(mydata,x='Treatment',y='logS_C',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Log(ARS)/Log(BG+XYL)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g10

#logS_N#
g11<-ggboxplot(mydata,x='Treatment',y='logS_N',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Log(ARS)/Log(NAG+LAP)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g11

#logS_P#
g12<-ggboxplot(mydata,x='Treatment',y='logS_P',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Log(ARS)/Log(AP)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g12

ggarrange(g9,g10,g11,g12,ncol = 2,nrow = 2,widths=c(1,1),heights=c(1,1),
common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)"))

#ggsave("plot2/Fig.3.pdf",plot = last_plot(),units = "in",height = 10,width = 14,dpi = 600)
####Fig. 4
#Fig. 4a
dataCON<-subset(mydata,Treatment=="CK")
dataCON.corr<-dataCON[,c(5:6,8:20)]
#controls of three forest
library(linkET)
qcorrplot(correlate(dataCON.corr,method = "spearman"), diag = FALSE,
grid_size = 0.25) +
geom_square() +
geom_mark(sep = '\n',size=5,sig_level = c(0.05,0.01,0.001),sig_thres = 0.05,
only_mark = TRUE)+
scale_fill_gradientn(colours = RColorBrewer::brewer.pal(11, "RdYlBu"),name="rho") +
guides(fill = guide_colorbar(title = "Spearman's rho", order = 3))

library(Hmisc)
res.cor<-rcorr(as.matrix(dataCON.corr),type = c("spearman"))
res.cor
## Total_C Total_N WSOC NH4 NO3 availP availS pH EC MBC MBN
## Total_C 1.00 0.82 0.77 0.48 0.78 0.74 0.45 -0.80 0.71 0.74 0.80
## Total_N 0.82 1.00 0.62 0.40 0.91 0.57 0.75 -0.87 0.85 0.92 0.89
## WSOC 0.77 0.62 1.00 0.60 0.66 0.60 0.25 -0.68 0.61 0.61 0.61
## NH4 0.48 0.40 0.60 1.00 0.33 0.20 0.25 -0.25 0.37 0.44 0.30
## NO3 0.78 0.91 0.66 0.33 1.00 0.71 0.56 -0.87 0.82 0.81 0.85
## availP 0.74 0.57 0.60 0.20 0.71 1.00 0.05 -0.62 0.53 0.41 0.50
## availS 0.45 0.75 0.25 0.25 0.56 0.05 1.00 -0.71 0.79 0.78 0.69
## pH -0.80 -0.87 -0.68 -0.25 -0.87 -0.62 -0.71 1.00 -0.96 -0.85 -0.87
## EC 0.71 0.85 0.61 0.37 0.82 0.53 0.79 -0.96 1.00 0.83 0.80
## MBC 0.74 0.92 0.61 0.44 0.81 0.41 0.78 -0.85 0.83 1.00 0.92
## MBN 0.80 0.89 0.61 0.30 0.85 0.50 0.69 -0.87 0.80 0.92 1.00
## MBC_MBN 0.30 0.45 0.45 0.63 0.34 0.10 0.49 -0.45 0.60 0.57 0.26
## ARS 0.70 0.62 0.66 0.79 0.63 0.51 0.25 -0.48 0.51 0.62 0.54
## SARSmbc -0.60 -0.81 -0.53 -0.22 -0.79 -0.37 -0.80 0.90 -0.87 -0.88 -0.89
## SARSt -0.67 -0.66 -0.64 0.01 -0.69 -0.58 -0.42 0.80 -0.65 -0.58 -0.68
## MBC_MBN ARS SARSmbc SARSt
## Total_C 0.30 0.70 -0.60 -0.67
## Total_N 0.45 0.62 -0.81 -0.66
## WSOC 0.45 0.66 -0.53 -0.64
## NH4 0.63 0.79 -0.22 0.01
## NO3 0.34 0.63 -0.79 -0.69
## availP 0.10 0.51 -0.37 -0.58
## availS 0.49 0.25 -0.80 -0.42
## pH -0.45 -0.48 0.90 0.80
## EC 0.60 0.51 -0.87 -0.65
## MBC 0.57 0.62 -0.88 -0.58
## MBN 0.26 0.54 -0.89 -0.68
## MBC_MBN 1.00 0.58 -0.39 -0.09
## ARS 0.58 1.00 -0.36 -0.11
## SARSmbc -0.39 -0.36 1.00 0.71
## SARSt -0.09 -0.11 0.71 1.00
##
## n= 15
##
##
## P
## Total_C Total_N WSOC NH4 NO3 availP availS pH EC MBC
## Total_C 0.0002 0.0008 0.0687 0.0006 0.0018 0.0953 0.0004 0.0028 0.0018
## Total_N 0.0002 0.0134 0.1358 0.0000 0.0261 0.0012 0.0000 0.0000 0.0000
## WSOC 0.0008 0.0134 0.0172 0.0069 0.0189 0.3760 0.0051 0.0148 0.0156
## NH4 0.0687 0.1358 0.0172 0.2265 0.4668 0.3618 0.3680 0.1728 0.1045
## NO3 0.0006 0.0000 0.0069 0.2265 0.0028 0.0297 0.0000 0.0002 0.0002
## availP 0.0018 0.0261 0.0189 0.4668 0.0028 0.8695 0.0142 0.0412 0.1283
## availS 0.0953 0.0012 0.3760 0.3618 0.0297 0.8695 0.0028 0.0005 0.0006
## pH 0.0004 0.0000 0.0051 0.3680 0.0000 0.0142 0.0028 0.0000 0.0000
## EC 0.0028 0.0000 0.0148 0.1728 0.0002 0.0412 0.0005 0.0000 0.0001
## MBC 0.0018 0.0000 0.0156 0.1045 0.0002 0.1283 0.0006 0.0000 0.0001
## MBN 0.0003 0.0000 0.0164 0.2773 0.0000 0.0598 0.0042 0.0000 0.0004 0.0000
## MBC_MBN 0.2834 0.0895 0.0924 0.0115 0.2212 0.7134 0.0664 0.0917 0.0189 0.0272
## ARS 0.0034 0.0127 0.0073 0.0005 0.0121 0.0537 0.3688 0.0718 0.0537 0.0134
## SARSmbc 0.0181 0.0003 0.0412 0.4277 0.0005 0.1773 0.0003 0.0000 0.0000 0.0000
## SARSt 0.0065 0.0073 0.0103 0.9597 0.0045 0.0238 0.1143 0.0004 0.0087 0.0238
## MBN MBC_MBN ARS SARSmbc SARSt
## Total_C 0.0003 0.2834 0.0034 0.0181 0.0065
## Total_N 0.0000 0.0895 0.0127 0.0003 0.0073
## WSOC 0.0164 0.0924 0.0073 0.0412 0.0103
## NH4 0.2773 0.0115 0.0005 0.4277 0.9597
## NO3 0.0000 0.2212 0.0121 0.0005 0.0045
## availP 0.0598 0.7134 0.0537 0.1773 0.0238
## availS 0.0042 0.0664 0.3688 0.0003 0.1143
## pH 0.0000 0.0917 0.0718 0.0000 0.0004
## EC 0.0004 0.0189 0.0537 0.0000 0.0087
## MBC 0.0000 0.0272 0.0134 0.0000 0.0238
## MBN 0.3412 0.0380 0.0000 0.0054
## MBC_MBN 0.3412 0.0228 0.1475 0.7613
## ARS 0.0380 0.0228 0.1866 0.6851
## SARSmbc 0.0000 0.1475 0.1866 0.0028
## SARSt 0.0054 0.7613 0.6851 0.0028
library(corrplot)#corrplot()#
#corrplot(res.cor$r,type = "upper",order = "original",p.mat = res.cor$P,sig.level = 0.05, insig = "blank",tl.pos = "tp")
#corrplot(res.cor$r,add = TRUE,type = "lower",method = "number",
#order = "original",number.cex = 0.6,diag = FALSE,tl.pos = "n",cl.pos = "n")
#Fig. 4b
dataCON<-mydata[mydata$Treatment=="CK",]
gr1<-ggplot(data = dataCON,aes(x=availS,y=ARS))+geom_point(aes(shape=Forest),size=5,color=c("#fdae61"))+
geom_smooth(method = "lm",se=FALSE)+
theme_bw()+ylab("Absolute ARS activity\n(nmol/g/h)")+xlab("Extractable S (mg/kg)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size),legend.position = "top")+
annotate("text",x=42,y=25,label="R2: 0.04\n p = 0.491",color="blue",size=5)
summary(lm(ARS~availS,data = dataCON))
##
## Call:
## lm(formula = ARS ~ availS, data = dataCON)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.921 -7.167 1.928 4.738 13.820
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.9227 9.1532 3.706 0.00264 **
## availS 0.2047 0.2890 0.708 0.49140
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.854 on 13 degrees of freedom
## Multiple R-squared: 0.03714, Adjusted R-squared: -0.03693
## F-statistic: 0.5014 on 1 and 13 DF, p-value: 0.4914
gr2<-ggplot(data = dataCON,aes(x=availS,y=SARSmbc))+geom_point(aes(shape=Forest),size=5,color=c("#fdae61"))+
geom_smooth(method = "lm",se=FALSE,color="red")+
theme_bw()+ylab("Specific ARS activity\n(nmol/ug MBC/h)")+xlab("Extractable S (mg/kg)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size),legend.position = "top")+
annotate("text",x=42,y=0.3,label="R2: 0.58\n p < 0.001",color="blue",size=5)
summary(lm(SARSmbc~availS,data = dataCON))
##
## Call:
## lm(formula = SARSmbc ~ availS, data = dataCON)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.09241 -0.03523 -0.01280 0.05062 0.09131
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.462208 0.057425 8.049 2.09e-06 ***
## availS -0.007714 0.001813 -4.254 0.00094 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06182 on 13 degrees of freedom
## Multiple R-squared: 0.5819, Adjusted R-squared: 0.5498
## F-statistic: 18.1 on 1 and 13 DF, p-value: 0.0009405
library(ggpubr)
ggarrange(gr1,gr2,ncol = 1,nrow = 2,widths=c(1,1),heights=c(1,1),
common.legend = TRUE,labels = c("(b)","(c)"))

#ggsave("plot2/Fig.4_availSenzyme.pdf",plot = last_plot(),units = "in",height = 7,width = 4,dpi = 600)
#Fig. 4c
####path analysis of dataCON####
library(lavaan)
library(semPlot)
dataCON.log<-mydata.log[mydata.log$Treatment == "CK",]
#ARS
#Fit good
model<-'
ARS~c1*MBC+c2*availS
MBC~b1*EC+b2*Total_C
availS~b3*availP+b4*EC
EC~a1*Total_C
availP~a2*Total_C
direct1:=a1*b1*c1
direct2:=a1*b4*c2
direct3:=a2*b3*c2
'
fit<-sem(model = model,data = dataCON.log,sample.nobs = 15)
summary(fit,standardized=TRUE,fit.measures=TRUE,rsq=TRUE)
## lavaan 0.6-12 ended normally after 2 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 13
##
## Number of observations 15
##
## Model Test User Model:
##
## Test statistic 8.826
## Degrees of freedom 7
## P-value (Chi-square) 0.265
##
## Model Test Baseline Model:
##
## Test statistic 109.611
## Degrees of freedom 15
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.981
## Tucker-Lewis Index (TLI) 0.959
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) 43.366
## Loglikelihood unrestricted model (H1) 47.779
##
## Akaike (AIC) -60.733
## Bayesian (BIC) -51.528
## Sample-size adjusted Bayesian (BIC) -91.216
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.132
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.361
## P-value RMSEA <= 0.05 0.290
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.068
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ARS ~
## MBC (c1) 0.538 0.100 5.383 0.000 0.538 1.058
## availS (c2) -0.526 0.173 -3.039 0.002 -0.526 -0.597
## MBC ~
## EC (b1) 1.140 0.214 5.334 0.000 1.140 0.813
## Total_C (b2) 0.218 0.204 1.071 0.284 0.218 0.163
## availS ~
## availP (b3) -1.463 0.484 -3.020 0.003 -1.463 -0.552
## EC (b4) 0.820 0.148 5.539 0.000 0.820 1.013
## EC ~
## Total_C (a1) 0.839 0.117 7.180 0.000 0.839 0.880
## availP ~
## Total_C (a2) 0.195 0.056 3.486 0.000 0.195 0.669
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .ARS 0.022 0.008 2.739 0.006 0.022 0.339
## .MBC 0.019 0.007 2.739 0.006 0.019 0.079
## .availS 0.027 0.010 2.739 0.006 0.027 0.328
## .EC 0.028 0.010 2.739 0.006 0.028 0.225
## .availP 0.007 0.002 2.739 0.006 0.007 0.552
##
## R-Square:
## Estimate
## ARS 0.661
## MBC 0.921
## availS 0.672
## EC 0.775
## availP 0.448
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## direct1 0.515 0.154 3.351 0.001 0.515 0.757
## direct2 -0.362 0.145 -2.498 0.012 -0.362 -0.532
## direct3 0.150 0.082 1.825 0.068 0.150 0.221
semPaths(fit,'std',layout = "tree2",sizeMan = 10,sizeLat = 10)

semPlot::semPaths(fit,'std',style="lisrel",layout = "tree2",
nodeLabels = c("ARS","MBC","Extractable S","EC","Extractable P","Total C"),
fade=FALSE,edge.label.cex = 1.75,sizeMan = 8,
residuals=FALSE,exoCov=FALSE)

#SARS
#Fit good
model<-'
SARSmbc~c1*MBC+c2*availS
MBC~b1*EC+b2*Total_C
availS~b3*availP+b4*EC
EC~a1*Total_C
availP~a2*Total_C
direct1:=a1*b1*c1
direct2:=a1*b4*c2
direct3:=a2*b3*c2
'
fit<-sem(model = model,data = dataCON.log,sample.nobs = 15)
summary(fit,standardized=TRUE,fit.measures=TRUE,rsq=TRUE)
## lavaan 0.6-12 ended normally after 3 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 13
##
## Number of observations 15
##
## Model Test User Model:
##
## Test statistic 8.930
## Degrees of freedom 7
## P-value (Chi-square) 0.258
##
## Model Test Baseline Model:
##
## Test statistic 119.494
## Degrees of freedom 15
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.982
## Tucker-Lewis Index (TLI) 0.960
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) 66.933
## Loglikelihood unrestricted model (H1) 71.398
##
## Akaike (AIC) -107.866
## Bayesian (BIC) -98.661
## Sample-size adjusted Bayesian (BIC) -138.349
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.136
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.363
## P-value RMSEA <= 0.05 0.282
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.037
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## SARSmbc ~
## MBC (c1) -0.073 0.021 -3.516 0.000 -0.073 -0.512
## availS (c2) -0.119 0.036 -3.323 0.001 -0.119 -0.484
## MBC ~
## EC (b1) 1.140 0.214 5.334 0.000 1.140 0.813
## Total_C (b2) 0.218 0.204 1.071 0.284 0.218 0.163
## availS ~
## availP (b3) -1.463 0.484 -3.020 0.003 -1.463 -0.552
## EC (b4) 0.820 0.148 5.539 0.000 0.820 1.013
## EC ~
## Total_C (a1) 0.839 0.117 7.180 0.000 0.839 0.880
## availP ~
## Total_C (a2) 0.195 0.056 3.486 0.000 0.195 0.669
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .SARSmbc 0.001 0.000 2.739 0.006 0.001 0.186
## .MBC 0.019 0.007 2.739 0.006 0.019 0.079
## .availS 0.027 0.010 2.739 0.006 0.027 0.328
## .EC 0.028 0.010 2.739 0.006 0.028 0.225
## .availP 0.007 0.002 2.739 0.006 0.007 0.552
##
## R-Square:
## Estimate
## SARSmbc 0.814
## MBC 0.921
## availS 0.672
## EC 0.775
## availP 0.448
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## direct1 -0.070 0.026 -2.717 0.007 -0.070 -0.366
## direct2 -0.082 0.031 -2.649 0.008 -0.082 -0.431
## direct3 0.034 0.018 1.881 0.060 0.034 0.179
semPaths(fit,'std',layout = "tree2",sizeMan = 10,sizeLat = 10)

semPlot::semPaths(fit,'std',style="lisrel",layout = "tree2",
nodeLabels = c("SARS","MBC","Extractable S","EC","Extractable P","Total C"),
fade=FALSE,edge.label.cex = 1.75,sizeMan = 8,
residuals=FALSE,exoCov=FALSE)

##Fig. 5
###Primary forest
library(psych)
library(tidyverse)
dataMF<-mydata[mydata$Forest == "PF",]
colnames(dataMF)
## [1] "Forest" "N" "P" "Treatment" "Total_C" "Total_N"
## [7] "TC_TN" "WSOC" "NH4" "NO3" "availP" "availS"
## [13] "pH" "EC" "MBC" "MBN" "MBC_MBN" "ARS"
## [19] "SARSmbc" "SARSt" "BG" "SBG" "NAG" "SNAG"
## [25] "XYL" "SXYL" "LAP" "SLAP" "ACP" "SACP"
## [31] "logS_C" "logS_N" "logS_P" "Forest2"
corr.tmp1<-dataMF[,c("Total_C","Total_N","TC_TN","WSOC","NH4","NO3","availP","availS",
"pH","EC","MBC","MBN","MBC_MBN")]
corr.tmp2<-dataMF[,c("ARS","SARSmbc","SARSt","logS_C","logS_N","logS_P")]
spman.d12 = corr.test(corr.tmp1, corr.tmp2,use="pairwise",method="spearman",adjust="fdr",alpha=.05,ci=FALSE)
r<-data.frame(spman.d12$r)
r[is.na(r)] <-0
r$X<-row.names(r)
r.long <- gather(r, Y, r, 1 : ncol(r)-1, factor_key=TRUE)
p<-data.frame(spman.d12$p)
p$X<-row.names(p)
p.long <- gather(p, Y, p, 1 : ncol(p)-1, factor_key=TRUE)
cor.out<-cbind(r.long,p.long$p)
cor.out$r <- round(as.numeric(cor.out$r), 2)
str(cor.out$r)
## num [1:78] -0.12 0.01 -0.09 0.39 0.07 -0.52 0.69 -0.69 0.26 0.13 ...
cor.out$X<-factor(cor.out$X,levels = c("Total_C","Total_N","TC_TN","WSOC","NH4","NO3","availP","availS",
"pH","EC","MBC","MBN","MBC_MBN"))
cor.out$X<-factor(cor.out$X,labels = c("Total C","Total N","Total C/total N","WSOC","NH4+","NO3-","extractable P","extractable S",
"pH","EC","MBC","MBN","MBC/MBN"))
cor.out$Y<-factor(cor.out$Y,levels = c("ARS","SARSmbc","SARSt","logS_C","logS_N","logS_P"))
cor.out$Y<-factor(cor.out$Y,labels = c("ARS","SARSmbc","SARStoc","LogS_C","LogS_N","LogS_P"))
ggplot(cor.out, aes(Y,X)) +
geom_tile(aes(fill = r), size=1)+
geom_text(aes(label = r), size = 1.5)+
scale_fill_gradient(guide = "legend", high='Salmon', low='CornflowerBlue')+
theme(axis.title= element_blank(),
axis.text.x=element_text(colour="black", size=8, face="bold",angle = 20),
axis.text.y=element_text(colour="black", size=8, face="bold"))

cor.out.sig <- cor.out
cor.out.sig$r[cor.out.sig$p > 0.05] <- 0
cor.out$Forest<-"PF"
cor.out.sig$Forest<-"PF"
ggplot(cor.out.sig, aes(Y,X)) +
geom_tile(aes(fill = r), size=1)+
geom_text(data= cor.out.sig[cor.out.sig$r!=0,], aes(label = r), size = 5, font="bold")+
scale_fill_gradient2(guide = "legend", high='#2c7bb6',mid="white", low='#d7191c')+
theme_bw()+
theme(axis.title= element_blank(),
axis.text.x=element_text(colour="black", size=15, face="bold",angle = 30, hjust=1),
axis.text.y=element_text(colour="black", size=15, face="bold"))

cor.out.PF<-cor.out
cor.out.sig.PF<-cor.out.sig
###Rehabilitated forest
library(psych)
library(tidyverse)
dataRF<-mydata[mydata$Forest=="RF",]
colnames(dataRF)
## [1] "Forest" "N" "P" "Treatment" "Total_C" "Total_N"
## [7] "TC_TN" "WSOC" "NH4" "NO3" "availP" "availS"
## [13] "pH" "EC" "MBC" "MBN" "MBC_MBN" "ARS"
## [19] "SARSmbc" "SARSt" "BG" "SBG" "NAG" "SNAG"
## [25] "XYL" "SXYL" "LAP" "SLAP" "ACP" "SACP"
## [31] "logS_C" "logS_N" "logS_P" "Forest2"
corr.tmp1<-dataRF[,c("Total_C","Total_N","TC_TN","WSOC","NH4","NO3","availP","availS",
"pH","EC","MBC","MBN","MBC_MBN")]
corr.tmp2<-dataRF[,c("ARS","SARSmbc","SARSt","logS_C","logS_N","logS_P")]
spman.d12 = corr.test(corr.tmp1, corr.tmp2,use="pairwise",method="spearman",adjust="fdr",alpha=.05,ci=FALSE)
r<-data.frame(spman.d12$r)
r[is.na(r)] <-0
r$X<-row.names(r)
r.long <- gather(r, Y, r, 1 : ncol(r)-1, factor_key=TRUE)
p<-data.frame(spman.d12$p)
p$X<-row.names(p)
p.long <- gather(p, Y, p, 1 : ncol(p)-1, factor_key=TRUE)
cor.out<-cbind(r.long,p.long$p)
cor.out$r <- round(as.numeric(cor.out$r), 2)
str(cor.out$r)
## num [1:78] 0.15 0.3 0.21 -0.18 -0.01 -0.01 0.08 -0.24 0.13 -0.15 ...
cor.out$X<-factor(cor.out$X,levels = c("Total_C","Total_N","TC_TN","WSOC","NH4","NO3","availP","availS",
"pH","EC","MBC","MBN","MBC_MBN"))
cor.out$X<-factor(cor.out$X,labels = c("Total C","Total N","Total C/total N","WSOC","NH4+","NO3-","extractable P","extractable S",
"pH","EC","MBC","MBN","MBC/MBN"))
cor.out$Y<-factor(cor.out$Y,levels = c("ARS","SARSmbc","SARSt","logS_C","logS_N","logS_P"))
cor.out$Y<-factor(cor.out$Y,labels = c("ARS","SARSmbc","SARStoc","LogS_C","LogS_N","LogS_P"))
ggplot(cor.out, aes(Y,X)) +
geom_tile(aes(fill = r), size=1)+
geom_text(aes(label = r), size = 1.5)+
scale_fill_gradient(guide = "legend", high='Salmon', low='CornflowerBlue')+
theme(axis.title= element_blank(),
axis.text.x=element_text(colour="black", size=8, face="bold",angle = 20),
axis.text.y=element_text(colour="black", size=8, face="bold"))

cor.out.sig <- cor.out
cor.out.sig$r[cor.out.sig$p > 0.05] <- 0
cor.out$Forest<-"RF"
cor.out.sig$Forest<-"RF"
ggplot(cor.out.sig, aes(Y,X)) +
geom_tile(aes(fill = r), size=1)+
geom_text(data= cor.out.sig[cor.out.sig$r!=0,], aes(label = r), size = 5, font="bold")+
scale_fill_gradient2(guide = "legend", high='#2c7bb6',mid="white", low='#d7191c')+
theme_bw()+
theme(axis.title= element_blank(),
axis.text.x=element_text(colour="black", size=15, face="bold",angle = 30, hjust=1),
axis.text.y=element_text(colour="black", size=15, face="bold"))

cor.out.RF<-cor.out
cor.out.sig.RF<-cor.out.sig
###Disturbed forest
library(psych)
library(tidyverse)
dataDF<-mydata[mydata$Forest == "DF",]
colnames(dataDF)
## [1] "Forest" "N" "P" "Treatment" "Total_C" "Total_N"
## [7] "TC_TN" "WSOC" "NH4" "NO3" "availP" "availS"
## [13] "pH" "EC" "MBC" "MBN" "MBC_MBN" "ARS"
## [19] "SARSmbc" "SARSt" "BG" "SBG" "NAG" "SNAG"
## [25] "XYL" "SXYL" "LAP" "SLAP" "ACP" "SACP"
## [31] "logS_C" "logS_N" "logS_P" "Forest2"
corr.tmp1<-dataDF[,c("Total_C","Total_N","TC_TN","WSOC","NH4","NO3","availP","availS",
"pH","EC","MBC","MBN","MBC_MBN")]
corr.tmp2<-dataDF[,c("ARS","SARSmbc","SARSt","logS_C","logS_N","logS_P")]
spman.d12 = corr.test(corr.tmp1, corr.tmp2,use="pairwise",method="spearman",adjust="fdr",alpha=.05,ci=FALSE)
r<-data.frame(spman.d12$r)
r[is.na(r)] <-0
r$X<-row.names(r)
r.long <- gather(r, Y, r, 1 : ncol(r)-1, factor_key=TRUE)
p<-data.frame(spman.d12$p)
p$X<-row.names(p)
p.long <- gather(p, Y, p, 1 : ncol(p)-1, factor_key=TRUE)
cor.out<-cbind(r.long,p.long$p)
cor.out$r <- round(as.numeric(cor.out$r), 2)
str(cor.out$r)
## num [1:78] 0.03 0.23 0.1 0.58 0.23 -0.41 0.72 -0.68 0.4 -0.25 ...
cor.out$X<-factor(cor.out$X,levels = c("Total_C","Total_N","TC_TN","WSOC","NH4","NO3","availP","availS",
"pH","EC","MBC","MBN","MBC_MBN"))
cor.out$X<-factor(cor.out$X,labels = c("Total C","Total N","Total C/total N","WSOC","NH4+","NO3-","extractable P","extractable S",
"pH","EC","MBC","MBN","MBC/MBN"))
cor.out$Y<-factor(cor.out$Y,levels = c("ARS","SARSmbc","SARSt","logS_C","logS_N","logS_P"))
cor.out$Y<-factor(cor.out$Y,labels = c("ARS","SARSmbc","SARStoc","LogS_C","LogS_N","LogS_P"))
ggplot(cor.out, aes(Y,X)) +
geom_tile(aes(fill = r), size=1)+
geom_text(aes(label = r), size = 1.5)+
scale_fill_gradient(guide = "legend", high='Salmon', low='CornflowerBlue')+
theme(axis.title= element_blank(),
axis.text.x=element_text(colour="black", size=8, face="bold",angle = 20),
axis.text.y=element_text(colour="black", size=8, face="bold"))

cor.out.sig <- cor.out
cor.out.sig$r[cor.out.sig$p > 0.05] <- 0
cor.out$Forest<-"DF"
cor.out.sig$Forest<-"DF"
ggplot(cor.out.sig, aes(Y,X)) +
geom_tile(aes(fill = r), size=1)+
geom_text(data= cor.out.sig[cor.out.sig$r!=0,], aes(label = r), size = 5, font="bold")+
scale_fill_gradient2(guide = "legend", high='#2c7bb6',mid="white", low='#d7191c')+
theme_bw()+
theme(axis.title= element_blank(),
axis.text.x=element_text(colour="black", size=15, face="bold",angle = 30, hjust=1),
axis.text.y=element_text(colour="black", size=15, face="bold"))

cor.out.DF<-cor.out
cor.out.sig.DF<-cor.out.sig
colnames(cor.out.PF) == colnames(cor.out.RF)
## [1] TRUE TRUE TRUE TRUE TRUE
colnames(cor.out.RF) == colnames(cor.out.DF)
## [1] TRUE TRUE TRUE TRUE TRUE
cor.out.all<-rbind(cor.out.PF,cor.out.RF,cor.out.DF)
cor.out.all$Forest<-factor(cor.out.all$Forest,levels = c("PF","RF","DF"),
labels = c("PF-Primary Forest","RF-Rehablitated Forest","DF-Disturbed Forest"))
colnames(cor.out.sig.PF) == colnames(cor.out.sig.RF)
## [1] TRUE TRUE TRUE TRUE TRUE
colnames(cor.out.sig.RF) == colnames(cor.out.sig.DF)
## [1] TRUE TRUE TRUE TRUE TRUE
cor.out.sig.all<-rbind(cor.out.sig.PF,cor.out.sig.RF,cor.out.sig.DF)
cor.out.sig.all$Forest<-factor(cor.out.sig.all$Forest,levels = c("PF","RF","DF"),
labels = c("PF-Primary Forest","RF-Rehablitated Forest","DF-Disturbed Forest"))
ggplot(cor.out.sig.all, aes(Y,X)) +
geom_tile(aes(fill = r), size=1)+
geom_text(data= cor.out.sig.all[cor.out.sig.all$r!=0,], aes(label = r), size = 5, font="bold")+
scale_fill_gradient2(guide = "legend", high='#2c7bb6',mid="white", low='#d7191c')+
theme_bw()+facet_grid(~Forest)+guides(fill = guide_legend(title = 'rho'))+
theme(axis.title= element_blank(),
axis.text.x=element_text(colour="black", size=15, face="bold",angle = 30, hjust=1),
axis.text.y=element_text(colour="black", size=15, face="bold"),
strip.text.x = element_text(size = 15,colour = "black"),
legend.key.size = unit(25,"pt"))

#ggsave("plot2/Fig.5_spearmanCorr.pdf",plot = last_plot(),units = "in",height = 6,width = 12,dpi = 600)
#Fig. 5b-c
####Data MF####
dataMF.log<-subset(mydata.log,Forest=="PF")
####stepwise regression analysis####
#ARS~
data.step<- dataMF.log[,5:18]
library(olsrr)
model<-ARS~.
fit<-lm(formula=model,data = data.step)
ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. Total_C
## 2. Total_N
## 3. TC_TN
## 4. WSOC
## 5. NH4
## 6. NO3
## 7. availP
## 8. availS
## 9. pH
## 10. EC
## 11. MBC
## 12. MBN
## 13. MBC_MBN
##
## We are selecting variables based on p value...
##
##
## Forward Selection: Step 1
##
## - availP
##
## Model Summary
## -------------------------------------------------------------
## R 0.820 RMSE 0.156
## R-Squared 0.672 Coef. Var 3.738
## Adj. R-Squared 0.654 MSE 0.024
## Pred R-Squared 0.592 MAE 0.110
## -------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 0.899 1 0.899 36.948 0.0000
## Residual 0.438 18 0.024
## Total 1.338 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 3.904 0.057 69.041 0.000 3.785 4.023
## availP 0.085 0.014 0.820 6.078 0.000 0.056 0.115
## -------------------------------------------------------------------------------------
##
##
##
## No more variables to be added.
##
## Variables Entered:
##
## + availP
##
##
## Final Model Output
## ------------------
##
## Model Summary
## -------------------------------------------------------------
## R 0.820 RMSE 0.156
## R-Squared 0.672 Coef. Var 3.738
## Adj. R-Squared 0.654 MSE 0.024
## Pred R-Squared 0.592 MAE 0.110
## -------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 0.899 1 0.899 36.948 0.0000
## Residual 0.438 18 0.024
## Total 1.338 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 3.904 0.057 69.041 0.000 3.785 4.023
## availP 0.085 0.014 0.820 6.078 0.000 0.056 0.115
## -------------------------------------------------------------------------------------
##
## Selection Summary
## -------------------------------------------------------------------------
## Variable Adj.
## Step Entered R-Square R-Square C(p) AIC RMSE
## -------------------------------------------------------------------------
## 1 availP 0.6724 0.6542 39.1403 -13.6600 0.1560
## -------------------------------------------------------------------------
test<-ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. Total_C
## 2. Total_N
## 3. TC_TN
## 4. WSOC
## 5. NH4
## 6. NO3
## 7. availP
## 8. availS
## 9. pH
## 10. EC
## 11. MBC
## 12. MBN
## 13. MBC_MBN
##
## We are selecting variables based on p value...
##
##
## Forward Selection: Step 1
##
## - availP
##
## Model Summary
## -------------------------------------------------------------
## R 0.820 RMSE 0.156
## R-Squared 0.672 Coef. Var 3.738
## Adj. R-Squared 0.654 MSE 0.024
## Pred R-Squared 0.592 MAE 0.110
## -------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 0.899 1 0.899 36.948 0.0000
## Residual 0.438 18 0.024
## Total 1.338 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 3.904 0.057 69.041 0.000 3.785 4.023
## availP 0.085 0.014 0.820 6.078 0.000 0.056 0.115
## -------------------------------------------------------------------------------------
##
##
##
## No more variables to be added.
##
## Variables Entered:
##
## + availP
##
##
## Final Model Output
## ------------------
##
## Model Summary
## -------------------------------------------------------------
## R 0.820 RMSE 0.156
## R-Squared 0.672 Coef. Var 3.738
## Adj. R-Squared 0.654 MSE 0.024
## Pred R-Squared 0.592 MAE 0.110
## -------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 0.899 1 0.899 36.948 0.0000
## Residual 0.438 18 0.024
## Total 1.338 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 3.904 0.057 69.041 0.000 3.785 4.023
## availP 0.085 0.014 0.820 6.078 0.000 0.056 0.115
## -------------------------------------------------------------------------------------
result.stepre.ars.MF<-data.frame(test$predictors,test$adjr)
colnames(result.stepre.ars.MF)<-c("prd","adjr")
#View(result.stepre.ars.MF)
result.stepre.ars.MF$prd2<-c("extractP")
result.stepre.ars.MF$adjp<-c("<0.001")
result.stepre.ars.MF$adjsig<-c("***")
result.stepre.ars.MF$adjr2<-c(0.6542)
library(ggplot2)
p1<-ggplot()+geom_col(data=result.stepre.ars.MF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
axis.title = element_text(size = 15,colour = "black"))+
geom_text(data = result.stepre.ars.MF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)
p1

#SARSm~
set.seed(123)
data.step<- dataMF.log[,c(5:17,19)]
library(olsrr)
model<-SARSmbc~.
fit<-lm(formula=model,data = data.step)
ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. Total_C
## 2. Total_N
## 3. TC_TN
## 4. WSOC
## 5. NH4
## 6. NO3
## 7. availP
## 8. availS
## 9. pH
## 10. EC
## 11. MBC
## 12. MBN
## 13. MBC_MBN
##
## We are selecting variables based on p value...
##
##
## Forward Selection: Step 1
##
## - MBC_MBN
##
## Model Summary
## --------------------------------------------------------------
## R 0.557 RMSE 0.031
## R-Squared 0.310 Coef. Var 18.784
## Adj. R-Squared 0.272 MSE 0.001
## Pred R-Squared 0.001 MAE 0.024
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.008 1 0.008 8.094 0.0107
## Residual 0.017 18 0.001
## Total 0.025 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 0.469 0.107 4.391 0.000 0.245 0.694
## MBC_MBN -0.166 0.058 -0.557 -2.845 0.011 -0.288 -0.043
## ----------------------------------------------------------------------------------------
##
##
##
## Forward Selection: Step 2
##
## - MBC
##
## Model Summary
## ---------------------------------------------------------------
## R 0.671 RMSE 0.029
## R-Squared 0.450 Coef. Var 17.260
## Adj. R-Squared 0.385 MSE 0.001
## Pred R-Squared -0.006 MAE 0.023
## ---------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.011 2 0.006 6.952 0.0062
## Residual 0.014 17 0.001
## Total 0.025 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 0.904 0.231 3.910 0.001 0.416 1.392
## MBC_MBN -0.183 0.054 -0.614 -3.373 0.004 -0.297 -0.068
## MBC -0.069 0.033 -0.378 -2.078 0.053 -0.138 0.001
## ----------------------------------------------------------------------------------------
##
##
##
## Forward Selection: Step 3
##
## - availP
##
## Model Summary
## --------------------------------------------------------------
## R 0.769 RMSE 0.025
## R-Squared 0.591 Coef. Var 15.345
## Adj. R-Squared 0.514 MSE 0.001
## Pred R-Squared 0.150 MAE 0.018
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.015 3 0.005 7.7 0.0021
## Residual 0.010 16 0.001
## Total 0.025 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 0.944 0.206 4.575 0.000 0.507 1.381
## MBC_MBN -0.043 0.077 -0.144 -0.559 0.584 -0.205 0.120
## MBC -0.124 0.038 -0.684 -3.293 0.005 -0.204 -0.044
## availP 0.010 0.004 0.706 2.347 0.032 0.001 0.019
## ----------------------------------------------------------------------------------------
##
##
##
## No more variables to be added.
##
## Variables Entered:
##
## + MBC_MBN
## + MBC
## + availP
##
##
## Final Model Output
## ------------------
##
## Model Summary
## --------------------------------------------------------------
## R 0.769 RMSE 0.025
## R-Squared 0.591 Coef. Var 15.345
## Adj. R-Squared 0.514 MSE 0.001
## Pred R-Squared 0.150 MAE 0.018
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.015 3 0.005 7.7 0.0021
## Residual 0.010 16 0.001
## Total 0.025 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 0.944 0.206 4.575 0.000 0.507 1.381
## MBC_MBN -0.043 0.077 -0.144 -0.559 0.584 -0.205 0.120
## MBC -0.124 0.038 -0.684 -3.293 0.005 -0.204 -0.044
## availP 0.010 0.004 0.706 2.347 0.032 0.001 0.019
## ----------------------------------------------------------------------------------------
##
## Selection Summary
## -------------------------------------------------------------------------
## Variable Adj.
## Step Entered R-Square R-Square C(p) AIC RMSE
## -------------------------------------------------------------------------
## 1 MBC_MBN 0.3102 0.2719 71.0093 -78.1093 0.0311
## 2 MBC 0.4499 0.3852 55.3836 -80.6366 0.0286
## 3 availP 0.5908 0.5141 39.6168 -84.5526 0.0254
## -------------------------------------------------------------------------
test<-ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. Total_C
## 2. Total_N
## 3. TC_TN
## 4. WSOC
## 5. NH4
## 6. NO3
## 7. availP
## 8. availS
## 9. pH
## 10. EC
## 11. MBC
## 12. MBN
## 13. MBC_MBN
##
## We are selecting variables based on p value...
##
##
## Forward Selection: Step 1
##
## - MBC_MBN
##
## Model Summary
## --------------------------------------------------------------
## R 0.557 RMSE 0.031
## R-Squared 0.310 Coef. Var 18.784
## Adj. R-Squared 0.272 MSE 0.001
## Pred R-Squared 0.001 MAE 0.024
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.008 1 0.008 8.094 0.0107
## Residual 0.017 18 0.001
## Total 0.025 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 0.469 0.107 4.391 0.000 0.245 0.694
## MBC_MBN -0.166 0.058 -0.557 -2.845 0.011 -0.288 -0.043
## ----------------------------------------------------------------------------------------
##
##
##
## Forward Selection: Step 2
##
## - MBC
##
## Model Summary
## ---------------------------------------------------------------
## R 0.671 RMSE 0.029
## R-Squared 0.450 Coef. Var 17.260
## Adj. R-Squared 0.385 MSE 0.001
## Pred R-Squared -0.006 MAE 0.023
## ---------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.011 2 0.006 6.952 0.0062
## Residual 0.014 17 0.001
## Total 0.025 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 0.904 0.231 3.910 0.001 0.416 1.392
## MBC_MBN -0.183 0.054 -0.614 -3.373 0.004 -0.297 -0.068
## MBC -0.069 0.033 -0.378 -2.078 0.053 -0.138 0.001
## ----------------------------------------------------------------------------------------
##
##
##
## Forward Selection: Step 3
##
## - availP
##
## Model Summary
## --------------------------------------------------------------
## R 0.769 RMSE 0.025
## R-Squared 0.591 Coef. Var 15.345
## Adj. R-Squared 0.514 MSE 0.001
## Pred R-Squared 0.150 MAE 0.018
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.015 3 0.005 7.7 0.0021
## Residual 0.010 16 0.001
## Total 0.025 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 0.944 0.206 4.575 0.000 0.507 1.381
## MBC_MBN -0.043 0.077 -0.144 -0.559 0.584 -0.205 0.120
## MBC -0.124 0.038 -0.684 -3.293 0.005 -0.204 -0.044
## availP 0.010 0.004 0.706 2.347 0.032 0.001 0.019
## ----------------------------------------------------------------------------------------
##
##
##
## No more variables to be added.
##
## Variables Entered:
##
## + MBC_MBN
## + MBC
## + availP
##
##
## Final Model Output
## ------------------
##
## Model Summary
## --------------------------------------------------------------
## R 0.769 RMSE 0.025
## R-Squared 0.591 Coef. Var 15.345
## Adj. R-Squared 0.514 MSE 0.001
## Pred R-Squared 0.150 MAE 0.018
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.015 3 0.005 7.7 0.0021
## Residual 0.010 16 0.001
## Total 0.025 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 0.944 0.206 4.575 0.000 0.507 1.381
## MBC_MBN -0.043 0.077 -0.144 -0.559 0.584 -0.205 0.120
## MBC -0.124 0.038 -0.684 -3.293 0.005 -0.204 -0.044
## availP 0.010 0.004 0.706 2.347 0.032 0.001 0.019
## ----------------------------------------------------------------------------------------
plot(test)


result.stepre.sars.MF<-data.frame(test$predictors,test$adjr)
colnames(result.stepre.sars.MF)<-c("prd","adjr")
#View(result.stepre.sars.MF)
result.stepre.sars.MF$prd2<-c("MBC_MBN","MBC","extractP")
result.stepre.sars.MF$adjp<-c("0.584","0.005","0.032")
result.stepre.sars.MF$adjsig<-c(".","***","*")
result.stepre.sars.MF$adjr2<-c(0.2719,0.1133,0.1289)
library(ggplot2)
p2<-ggplot()+geom_col(data=result.stepre.sars.MF,aes(x=prd2,y=adjr2),fill="#fdae61",color=NA,width = 0.5)+
labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
axis.title = element_text(size = 15,colour = "black"))+
geom_text(data = result.stepre.sars.MF,aes(x=prd2,y=adjr2,label=adjsig),nudge_y = 0.05,size=7)
p2

library(ggpubr)
ggarrange(p1,p2,ncol = 1,nrow = 2,widths=c(1,1),heights=c(1,1),labels = c("(b)","(c)"))

#ggsave("plot2/Fig.2_stepwiseControls.pdf",plot = last_plot(),units = "in",height = 12,width = 6,dpi = 600)
result.stepre.ars.MF$enzyme<-"ARS-PF"
result.stepre.sars.MF$enzyme<-"SARS-PF"
result.stepre.MF<-rbind(result.stepre.ars.MF,result.stepre.sars.MF)
result.stepre.MF$forest<-"PF"
#View(result.stepre.MF)
library(ggplot2)
p3<-ggplot()+geom_col(data=result.stepre.MF,aes(x=prd2,y=adjr2),fill="#fdae61",color=NA,width = 0.5)+
labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
axis.title = element_text(size = 15,colour = "black"),
axis.text.x = element_text(angle = 15,hjust = 1),
strip.text.x = element_text(size = 15,colour = "black"))+
geom_text(data = result.stepre.MF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)+
facet_grid(~enzyme,scales = "free")
p3

####Data RF####
dataRF.log<-subset(mydata.log,Forest=="RF")
####stepwise regression analysis####
#ARS~
data.step<- dataRF.log[,5:18]
library(olsrr)
model<-ARS~.
fit<-lm(formula=model,data = data.step)
ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. Total_C
## 2. Total_N
## 3. TC_TN
## 4. WSOC
## 5. NH4
## 6. NO3
## 7. availP
## 8. availS
## 9. pH
## 10. EC
## 11. MBC
## 12. MBN
## 13. MBC_MBN
##
## We are selecting variables based on p value...
##
##
## Forward Selection: Step 1
##
## - MBC
##
## Model Summary
## --------------------------------------------------------------
## R 0.378 RMSE 0.168
## R-Squared 0.143 Coef. Var 4.617
## Adj. R-Squared 0.096 MSE 0.028
## Pred R-Squared -0.024 MAE 0.139
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.085 1 0.085 3.009 0.0999
## Residual 0.507 18 0.028
## Total 0.592 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 2.349 0.743 3.162 0.005 0.788 3.910
## MBC 0.264 0.152 0.378 1.735 0.100 -0.056 0.584
## -------------------------------------------------------------------------------------
##
##
##
## No more variables to be added.
##
## Variables Entered:
##
## + MBC
##
##
## Final Model Output
## ------------------
##
## Model Summary
## --------------------------------------------------------------
## R 0.378 RMSE 0.168
## R-Squared 0.143 Coef. Var 4.617
## Adj. R-Squared 0.096 MSE 0.028
## Pred R-Squared -0.024 MAE 0.139
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.085 1 0.085 3.009 0.0999
## Residual 0.507 18 0.028
## Total 0.592 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 2.349 0.743 3.162 0.005 0.788 3.910
## MBC 0.264 0.152 0.378 1.735 0.100 -0.056 0.584
## -------------------------------------------------------------------------------------
##
## Selection Summary
## ------------------------------------------------------------------------
## Variable Adj.
## Step Entered R-Square R-Square C(p) AIC RMSE
## ------------------------------------------------------------------------
## 1 MBC 0.1432 0.0956 9.1247 -10.7299 0.1679
## ------------------------------------------------------------------------
test<-ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. Total_C
## 2. Total_N
## 3. TC_TN
## 4. WSOC
## 5. NH4
## 6. NO3
## 7. availP
## 8. availS
## 9. pH
## 10. EC
## 11. MBC
## 12. MBN
## 13. MBC_MBN
##
## We are selecting variables based on p value...
##
##
## Forward Selection: Step 1
##
## - MBC
##
## Model Summary
## --------------------------------------------------------------
## R 0.378 RMSE 0.168
## R-Squared 0.143 Coef. Var 4.617
## Adj. R-Squared 0.096 MSE 0.028
## Pred R-Squared -0.024 MAE 0.139
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.085 1 0.085 3.009 0.0999
## Residual 0.507 18 0.028
## Total 0.592 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 2.349 0.743 3.162 0.005 0.788 3.910
## MBC 0.264 0.152 0.378 1.735 0.100 -0.056 0.584
## -------------------------------------------------------------------------------------
##
##
##
## No more variables to be added.
##
## Variables Entered:
##
## + MBC
##
##
## Final Model Output
## ------------------
##
## Model Summary
## --------------------------------------------------------------
## R 0.378 RMSE 0.168
## R-Squared 0.143 Coef. Var 4.617
## Adj. R-Squared 0.096 MSE 0.028
## Pred R-Squared -0.024 MAE 0.139
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.085 1 0.085 3.009 0.0999
## Residual 0.507 18 0.028
## Total 0.592 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 2.349 0.743 3.162 0.005 0.788 3.910
## MBC 0.264 0.152 0.378 1.735 0.100 -0.056 0.584
## -------------------------------------------------------------------------------------
result.stepre.ars.RF<-data.frame(test$predictors,test$adjr)
colnames(result.stepre.ars.RF)<-c("prd","adjr")
#View(result.stepre.ars.RF)
result.stepre.ars.RF$prd2<-c("MBC")
result.stepre.ars.RF$adjp<-c("0.100")
result.stepre.ars.RF$adjsig<-c(".")
result.stepre.ars.RF$adjr2<-c(0.1432)
library(ggplot2)
p4<-ggplot()+geom_col(data=result.stepre.ars.RF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
axis.title = element_text(size = 15,colour = "black"))+
geom_text(data = result.stepre.ars.RF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)
p4

#SARSm~
data.step<- dataRF.log[,c(5:17,19)]#dataCON[,c(3:15,17)]
library(olsrr)
model<-SARSmbc~.
fit<-lm(formula=model,data = data.step)
ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. Total_C
## 2. Total_N
## 3. TC_TN
## 4. WSOC
## 5. NH4
## 6. NO3
## 7. availP
## 8. availS
## 9. pH
## 10. EC
## 11. MBC
## 12. MBN
## 13. MBC_MBN
##
## We are selecting variables based on p value...
##
##
## Forward Selection: Step 1
##
## - MBC
##
## Model Summary
## --------------------------------------------------------------
## R 0.748 RMSE 0.039
## R-Squared 0.560 Coef. Var 15.306
## Adj. R-Squared 0.536 MSE 0.002
## Pred R-Squared 0.463 MAE 0.033
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 0.035 1 0.035 22.916 1e-04
## Residual 0.028 18 0.002
## Total 0.063 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 1.083 0.173 6.257 0.000 0.719 1.446
## MBC -0.170 0.035 -0.748 -4.787 0.000 -0.244 -0.095
## ----------------------------------------------------------------------------------------
##
##
##
## No more variables to be added.
##
## Variables Entered:
##
## + MBC
##
##
## Final Model Output
## ------------------
##
## Model Summary
## --------------------------------------------------------------
## R 0.748 RMSE 0.039
## R-Squared 0.560 Coef. Var 15.306
## Adj. R-Squared 0.536 MSE 0.002
## Pred R-Squared 0.463 MAE 0.033
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 0.035 1 0.035 22.916 1e-04
## Residual 0.028 18 0.002
## Total 0.063 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 1.083 0.173 6.257 0.000 0.719 1.446
## MBC -0.170 0.035 -0.748 -4.787 0.000 -0.244 -0.095
## ----------------------------------------------------------------------------------------
##
## Selection Summary
## ------------------------------------------------------------------------
## Variable Adj.
## Step Entered R-Square R-Square C(p) AIC RMSE
## ------------------------------------------------------------------------
## 1 MBC 0.5601 0.5356 7.0863 -69.0077 0.0391
## ------------------------------------------------------------------------
test<-ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. Total_C
## 2. Total_N
## 3. TC_TN
## 4. WSOC
## 5. NH4
## 6. NO3
## 7. availP
## 8. availS
## 9. pH
## 10. EC
## 11. MBC
## 12. MBN
## 13. MBC_MBN
##
## We are selecting variables based on p value...
##
##
## Forward Selection: Step 1
##
## - MBC
##
## Model Summary
## --------------------------------------------------------------
## R 0.748 RMSE 0.039
## R-Squared 0.560 Coef. Var 15.306
## Adj. R-Squared 0.536 MSE 0.002
## Pred R-Squared 0.463 MAE 0.033
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 0.035 1 0.035 22.916 1e-04
## Residual 0.028 18 0.002
## Total 0.063 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 1.083 0.173 6.257 0.000 0.719 1.446
## MBC -0.170 0.035 -0.748 -4.787 0.000 -0.244 -0.095
## ----------------------------------------------------------------------------------------
##
##
##
## No more variables to be added.
##
## Variables Entered:
##
## + MBC
##
##
## Final Model Output
## ------------------
##
## Model Summary
## --------------------------------------------------------------
## R 0.748 RMSE 0.039
## R-Squared 0.560 Coef. Var 15.306
## Adj. R-Squared 0.536 MSE 0.002
## Pred R-Squared 0.463 MAE 0.033
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 0.035 1 0.035 22.916 1e-04
## Residual 0.028 18 0.002
## Total 0.063 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 1.083 0.173 6.257 0.000 0.719 1.446
## MBC -0.170 0.035 -0.748 -4.787 0.000 -0.244 -0.095
## ----------------------------------------------------------------------------------------
result.stepre.sars.RF<-data.frame(test$predictors,test$adjr)
colnames(result.stepre.sars.RF)<-c("prd","adjr")
#View(result.stepre.sars.RF)
result.stepre.sars.RF$prd2<-c("MBC")
result.stepre.sars.RF$adjp<-c("0.0001")
result.stepre.sars.RF$adjsig<-c("***")
result.stepre.sars.RF$adjr2<-c(0.5356)
library(ggplot2)
p5<-ggplot()+geom_col(data=result.stepre.sars.RF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
axis.title = element_text(size = 15,colour = "black"))+
geom_text(data = result.stepre.sars.RF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)
p5

library(ggpubr)
ggarrange(p1,p2,ncol = 1,nrow = 2,widths=c(1,1),heights=c(1,1),labels = c("(b)","(c)"))

#ggsave("plot2/Fig.2_stepwiseControls.pdf",plot = last_plot(),units = "in",height = 12,width = 6,dpi = 600)
result.stepre.ars.RF$enzyme<-"ARS-RF"
result.stepre.sars.RF$enzyme<-"SARS-RF"
result.stepre.RF<-rbind(result.stepre.ars.RF,result.stepre.sars.RF)
result.stepre.RF$forest<-"RF"
#View(result.stepre.RF)
library(ggplot2)
p6<-ggplot()+geom_col(data=result.stepre.RF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
axis.title = element_text(size = 15,colour = "black"),
#axis.text.x = element_text(angle = 15,hjust = 1),
strip.text.x = element_text(size = 15,colour = "black"))+
geom_text(data = result.stepre.RF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)+
facet_grid(~enzyme,scales = "free")
p6

####Data DF####
dataDF.log<-subset(mydata.log,Forest=="DF")
####stepwise regression analysis####
#ARS~
data.step<- dataDF.log[,5:18]
library(olsrr)
model<-ARS~.
fit<-lm(formula=model,data = data.step)
ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. Total_C
## 2. Total_N
## 3. TC_TN
## 4. WSOC
## 5. NH4
## 6. NO3
## 7. availP
## 8. availS
## 9. pH
## 10. EC
## 11. MBC
## 12. MBN
## 13. MBC_MBN
##
## We are selecting variables based on p value...
##
##
## Forward Selection: Step 1
##
## - availP
##
## Model Summary
## -------------------------------------------------------------
## R 0.811 RMSE 0.230
## R-Squared 0.658 Coef. Var 6.121
## Adj. R-Squared 0.639 MSE 0.053
## Pred R-Squared 0.590 MAE 0.168
## -------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 1.837 1 1.837 34.581 0.0000
## Residual 0.956 18 0.053
## Total 2.793 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 3.385 0.083 40.902 0.000 3.211 3.559
## availP 0.123 0.021 0.811 5.881 0.000 0.079 0.166
## -------------------------------------------------------------------------------------
##
##
##
## Forward Selection: Step 2
##
## - NH4
##
## Model Summary
## -------------------------------------------------------------
## R 0.850 RMSE 0.214
## R-Squared 0.722 Coef. Var 5.676
## Adj. R-Squared 0.689 MSE 0.046
## Pred R-Squared 0.640 MAE 0.162
## -------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 2.017 2 1.008 22.079 0.0000
## Residual 0.776 17 0.046
## Total 2.793 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 2.773 0.318 8.721 0.000 2.102 3.443
## availP 0.121 0.019 0.802 6.268 0.000 0.080 0.162
## NH4 0.726 0.366 0.254 1.984 0.064 -0.046 1.498
## -------------------------------------------------------------------------------------
##
##
##
## No more variables to be added.
##
## Variables Entered:
##
## + availP
## + NH4
##
##
## Final Model Output
## ------------------
##
## Model Summary
## -------------------------------------------------------------
## R 0.850 RMSE 0.214
## R-Squared 0.722 Coef. Var 5.676
## Adj. R-Squared 0.689 MSE 0.046
## Pred R-Squared 0.640 MAE 0.162
## -------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 2.017 2 1.008 22.079 0.0000
## Residual 0.776 17 0.046
## Total 2.793 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 2.773 0.318 8.721 0.000 2.102 3.443
## availP 0.121 0.019 0.802 6.268 0.000 0.080 0.162
## NH4 0.726 0.366 0.254 1.984 0.064 -0.046 1.498
## -------------------------------------------------------------------------------------
##
## Selection Summary
## ------------------------------------------------------------------------
## Variable Adj.
## Step Entered R-Square R-Square C(p) AIC RMSE
## ------------------------------------------------------------------------
## 1 availP 0.6577 0.6387 -2.5829 1.9478 0.2305
## 2 NH4 0.7220 0.6893 -3.1055 -0.2176 0.2137
## ------------------------------------------------------------------------
test<-ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. Total_C
## 2. Total_N
## 3. TC_TN
## 4. WSOC
## 5. NH4
## 6. NO3
## 7. availP
## 8. availS
## 9. pH
## 10. EC
## 11. MBC
## 12. MBN
## 13. MBC_MBN
##
## We are selecting variables based on p value...
##
##
## Forward Selection: Step 1
##
## - availP
##
## Model Summary
## -------------------------------------------------------------
## R 0.811 RMSE 0.230
## R-Squared 0.658 Coef. Var 6.121
## Adj. R-Squared 0.639 MSE 0.053
## Pred R-Squared 0.590 MAE 0.168
## -------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 1.837 1 1.837 34.581 0.0000
## Residual 0.956 18 0.053
## Total 2.793 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 3.385 0.083 40.902 0.000 3.211 3.559
## availP 0.123 0.021 0.811 5.881 0.000 0.079 0.166
## -------------------------------------------------------------------------------------
##
##
##
## Forward Selection: Step 2
##
## - NH4
##
## Model Summary
## -------------------------------------------------------------
## R 0.850 RMSE 0.214
## R-Squared 0.722 Coef. Var 5.676
## Adj. R-Squared 0.689 MSE 0.046
## Pred R-Squared 0.640 MAE 0.162
## -------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 2.017 2 1.008 22.079 0.0000
## Residual 0.776 17 0.046
## Total 2.793 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 2.773 0.318 8.721 0.000 2.102 3.443
## availP 0.121 0.019 0.802 6.268 0.000 0.080 0.162
## NH4 0.726 0.366 0.254 1.984 0.064 -0.046 1.498
## -------------------------------------------------------------------------------------
##
##
##
## No more variables to be added.
##
## Variables Entered:
##
## + availP
## + NH4
##
##
## Final Model Output
## ------------------
##
## Model Summary
## -------------------------------------------------------------
## R 0.850 RMSE 0.214
## R-Squared 0.722 Coef. Var 5.676
## Adj. R-Squared 0.689 MSE 0.046
## Pred R-Squared 0.640 MAE 0.162
## -------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 2.017 2 1.008 22.079 0.0000
## Residual 0.776 17 0.046
## Total 2.793 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------
## (Intercept) 2.773 0.318 8.721 0.000 2.102 3.443
## availP 0.121 0.019 0.802 6.268 0.000 0.080 0.162
## NH4 0.726 0.366 0.254 1.984 0.064 -0.046 1.498
## -------------------------------------------------------------------------------------
result.stepre.ars.DF<-data.frame(test$predictors,test$adjr)
colnames(result.stepre.ars.DF)<-c("prd","adjr")
#View(result.stepre.ars.DF)
result.stepre.ars.DF$prd2<-c("extactP","NH4+")
result.stepre.ars.DF$adjp<-c("0.0001","0.064")
result.stepre.ars.DF$adjsig<-c("***",".")
result.stepre.ars.DF$adjr2<-c(0.6387,0.0506)
library(ggplot2)
p7<-ggplot()+geom_col(data=result.stepre.ars.DF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
axis.title = element_text(size = 15,colour = "black"))+
geom_text(data = result.stepre.ars.DF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)
p7

#SARSm~
data.step<- dataDF.log[,c(5:17,19)]#dataCON[,c(3:15,17)]
library(olsrr)
model<-SARSmbc~.
fit<-lm(formula=model,data = data.step)
ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. Total_C
## 2. Total_N
## 3. TC_TN
## 4. WSOC
## 5. NH4
## 6. NO3
## 7. availP
## 8. availS
## 9. pH
## 10. EC
## 11. MBC
## 12. MBN
## 13. MBC_MBN
##
## We are selecting variables based on p value...
##
##
## Forward Selection: Step 1
##
## - availS
##
## Model Summary
## --------------------------------------------------------------
## R 0.594 RMSE 0.058
## R-Squared 0.352 Coef. Var 25.110
## Adj. R-Squared 0.316 MSE 0.003
## Pred R-Squared 0.210 MAE 0.047
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.033 1 0.033 9.793 0.0058
## Residual 0.061 18 0.003
## Total 0.095 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 0.664 0.139 4.794 0.000 0.373 0.955
## availS -0.131 0.042 -0.594 -3.129 0.006 -0.219 -0.043
## ----------------------------------------------------------------------------------------
##
##
##
## Forward Selection: Step 2
##
## - MBC
##
## Model Summary
## --------------------------------------------------------------
## R 0.694 RMSE 0.054
## R-Squared 0.482 Coef. Var 23.103
## Adj. R-Squared 0.421 MSE 0.003
## Pred R-Squared 0.271 MAE 0.041
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.046 2 0.023 7.916 0.0037
## Residual 0.049 17 0.003
## Total 0.095 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 1.403 0.380 3.694 0.002 0.602 2.204
## availS -0.153 0.040 -0.695 -3.834 0.001 -0.238 -0.069
## MBC -0.129 0.063 -0.374 -2.065 0.055 -0.262 0.003
## ----------------------------------------------------------------------------------------
##
##
##
## Forward Selection: Step 3
##
## - availP
##
## Model Summary
## --------------------------------------------------------------
## R 0.816 RMSE 0.044
## R-Squared 0.666 Coef. Var 19.134
## Adj. R-Squared 0.603 MSE 0.002
## Pred R-Squared 0.454 MAE 0.031
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 0.063 3 0.021 10.621 4e-04
## Residual 0.032 16 0.002
## Total 0.095 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 1.465 0.315 4.649 0.000 0.797 2.133
## availS -0.046 0.049 -0.207 -0.929 0.367 -0.150 0.059
## MBC -0.223 0.061 -0.646 -3.673 0.002 -0.352 -0.094
## availP 0.021 0.007 0.754 2.964 0.009 0.006 0.036
## ----------------------------------------------------------------------------------------
##
##
##
## Forward Selection: Step 4
##
## - NH4
##
## Model Summary
## --------------------------------------------------------------
## R 0.851 RMSE 0.042
## R-Squared 0.724 Coef. Var 17.946
## Adj. R-Squared 0.651 MSE 0.002
## Pred R-Squared 0.499 MAE 0.030
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.069 4 0.017 9.853 4e-04
## Residual 0.026 15 0.002
## Total 0.095 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 1.343 0.303 4.429 0.000 0.697 1.990
## availS -0.018 0.049 -0.081 -0.369 0.717 -0.122 0.086
## MBC -0.242 0.058 -0.700 -4.171 0.001 -0.365 -0.118
## availP 0.024 0.007 0.867 3.513 0.003 0.009 0.039
## NH4 0.135 0.076 0.257 1.786 0.094 -0.026 0.297
## ----------------------------------------------------------------------------------------
##
##
##
## No more variables to be added.
##
## Variables Entered:
##
## + availS
## + MBC
## + availP
## + NH4
##
##
## Final Model Output
## ------------------
##
## Model Summary
## --------------------------------------------------------------
## R 0.851 RMSE 0.042
## R-Squared 0.724 Coef. Var 17.946
## Adj. R-Squared 0.651 MSE 0.002
## Pred R-Squared 0.499 MAE 0.030
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.069 4 0.017 9.853 4e-04
## Residual 0.026 15 0.002
## Total 0.095 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 1.343 0.303 4.429 0.000 0.697 1.990
## availS -0.018 0.049 -0.081 -0.369 0.717 -0.122 0.086
## MBC -0.242 0.058 -0.700 -4.171 0.001 -0.365 -0.118
## availP 0.024 0.007 0.867 3.513 0.003 0.009 0.039
## NH4 0.135 0.076 0.257 1.786 0.094 -0.026 0.297
## ----------------------------------------------------------------------------------------
##
## Selection Summary
## ------------------------------------------------------------------------
## Variable Adj.
## Step Entered R-Square R-Square C(p) AIC RMSE
## ------------------------------------------------------------------------
## 1 availS 0.3523 0.3164 8.7795 -52.9700 0.0584
## 2 MBC 0.4822 0.4213 5.8109 -55.4457 0.0537
## 3 availP 0.6657 0.6030 0.7898 -62.1974 0.0445
## 4 NH4 0.7243 0.6508 0.5477 -64.0521 0.0417
## ------------------------------------------------------------------------
test<-ols_step_forward_p(fit, pent=0.1,perm=0.3,details = TRUE)
## Forward Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. Total_C
## 2. Total_N
## 3. TC_TN
## 4. WSOC
## 5. NH4
## 6. NO3
## 7. availP
## 8. availS
## 9. pH
## 10. EC
## 11. MBC
## 12. MBN
## 13. MBC_MBN
##
## We are selecting variables based on p value...
##
##
## Forward Selection: Step 1
##
## - availS
##
## Model Summary
## --------------------------------------------------------------
## R 0.594 RMSE 0.058
## R-Squared 0.352 Coef. Var 25.110
## Adj. R-Squared 0.316 MSE 0.003
## Pred R-Squared 0.210 MAE 0.047
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.033 1 0.033 9.793 0.0058
## Residual 0.061 18 0.003
## Total 0.095 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 0.664 0.139 4.794 0.000 0.373 0.955
## availS -0.131 0.042 -0.594 -3.129 0.006 -0.219 -0.043
## ----------------------------------------------------------------------------------------
##
##
##
## Forward Selection: Step 2
##
## - MBC
##
## Model Summary
## --------------------------------------------------------------
## R 0.694 RMSE 0.054
## R-Squared 0.482 Coef. Var 23.103
## Adj. R-Squared 0.421 MSE 0.003
## Pred R-Squared 0.271 MAE 0.041
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.046 2 0.023 7.916 0.0037
## Residual 0.049 17 0.003
## Total 0.095 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 1.403 0.380 3.694 0.002 0.602 2.204
## availS -0.153 0.040 -0.695 -3.834 0.001 -0.238 -0.069
## MBC -0.129 0.063 -0.374 -2.065 0.055 -0.262 0.003
## ----------------------------------------------------------------------------------------
##
##
##
## Forward Selection: Step 3
##
## - availP
##
## Model Summary
## --------------------------------------------------------------
## R 0.816 RMSE 0.044
## R-Squared 0.666 Coef. Var 19.134
## Adj. R-Squared 0.603 MSE 0.002
## Pred R-Squared 0.454 MAE 0.031
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## -------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -------------------------------------------------------------------
## Regression 0.063 3 0.021 10.621 4e-04
## Residual 0.032 16 0.002
## Total 0.095 19
## -------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 1.465 0.315 4.649 0.000 0.797 2.133
## availS -0.046 0.049 -0.207 -0.929 0.367 -0.150 0.059
## MBC -0.223 0.061 -0.646 -3.673 0.002 -0.352 -0.094
## availP 0.021 0.007 0.754 2.964 0.009 0.006 0.036
## ----------------------------------------------------------------------------------------
##
##
##
## Forward Selection: Step 4
##
## - NH4
##
## Model Summary
## --------------------------------------------------------------
## R 0.851 RMSE 0.042
## R-Squared 0.724 Coef. Var 17.946
## Adj. R-Squared 0.651 MSE 0.002
## Pred R-Squared 0.499 MAE 0.030
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.069 4 0.017 9.853 4e-04
## Residual 0.026 15 0.002
## Total 0.095 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 1.343 0.303 4.429 0.000 0.697 1.990
## availS -0.018 0.049 -0.081 -0.369 0.717 -0.122 0.086
## MBC -0.242 0.058 -0.700 -4.171 0.001 -0.365 -0.118
## availP 0.024 0.007 0.867 3.513 0.003 0.009 0.039
## NH4 0.135 0.076 0.257 1.786 0.094 -0.026 0.297
## ----------------------------------------------------------------------------------------
##
##
##
## No more variables to be added.
##
## Variables Entered:
##
## + availS
## + MBC
## + availP
## + NH4
##
##
## Final Model Output
## ------------------
##
## Model Summary
## --------------------------------------------------------------
## R 0.851 RMSE 0.042
## R-Squared 0.724 Coef. Var 17.946
## Adj. R-Squared 0.651 MSE 0.002
## Pred R-Squared 0.499 MAE 0.030
## --------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.069 4 0.017 9.853 4e-04
## Residual 0.026 15 0.002
## Total 0.095 19
## ------------------------------------------------------------------
##
## Parameter Estimates
## ----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## ----------------------------------------------------------------------------------------
## (Intercept) 1.343 0.303 4.429 0.000 0.697 1.990
## availS -0.018 0.049 -0.081 -0.369 0.717 -0.122 0.086
## MBC -0.242 0.058 -0.700 -4.171 0.001 -0.365 -0.118
## availP 0.024 0.007 0.867 3.513 0.003 0.009 0.039
## NH4 0.135 0.076 0.257 1.786 0.094 -0.026 0.297
## ----------------------------------------------------------------------------------------
result.stepre.sars.DF<-data.frame(test$predictors,test$adjr)
colnames(result.stepre.sars.DF)<-c("prd","adjr")
#View(result.stepre.sars.DF)
result.stepre.sars.DF$prd2<-c("extractS","MBC","extractP","NH4")
result.stepre.sars.DF$adjp<-c("0.717","0.001","0.003","0.094")
result.stepre.sars.DF$adjsig<-c(".","***","**",".")
result.stepre.sars.DF$adjr2<-c(0.3164,0.1049,0.1817,0.0478)
result.stepre.sars.DF<-result.stepre.sars.DF[result.stepre.sars.DF$adjp<0.05,]
library(ggplot2)
p8<-ggplot()+geom_col(data=result.stepre.sars.DF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
axis.title = element_text(size = 15,colour = "black"))+
geom_text(data = result.stepre.sars.DF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)
p8

library(ggpubr)
ggarrange(p1,p2,ncol = 1,nrow = 2,widths=c(1,1),heights=c(1,1),labels = c("(b)","(c)"))

#ggsave("plot2/Fig.2_stepwiseControls.pdf",plot = last_plot(),units = "in",height = 12,width = 6,dpi = 600)
result.stepre.ars.DF$enzyme<-"ARS-DF"
result.stepre.sars.DF$enzyme<-"SARS-DF"
result.stepre.DF<-rbind(result.stepre.ars.DF,result.stepre.sars.DF)
result.stepre.DF$forest<-"DF"
#View(result.stepre.DF)
library(ggplot2)
p9<-ggplot()+geom_col(data=result.stepre.DF,aes(x=prd2,y=adjr),fill="#fdae61",color=NA,width = 0.5)+
labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
axis.title = element_text(size = 15,colour = "black"),
#axis.text.x = element_text(angle = 15,hjust =1),
strip.text.x = element_text(size = 15,colour = "black"))+
geom_text(data = result.stepre.DF,aes(x=prd2,y=adjr,label=adjsig),nudge_y = 0.05,size=7)+
facet_grid(~enzyme,scales = "free")
p9

library(ggpubr)
ggarrange(p3,p6,p9,ncol = 3,nrow = 1,widths=c(1,1),heights=c(1,1),labels = c("(b)","(c)","(d)"))

#ggsave("plot2/Fig.5_stepwiseControls.pdf",plot = last_plot(),units = "in",height = 4,width = 12,dpi = 600)
library(ggplot2)
p3<-ggplot(data=result.stepre.MF,aes(x=reorder(prd2, -adjr2),y=adjr2))+geom_col(fill="#fdae61",color=NA,width = 0.5)+
labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
axis.title = element_text(size = 15,colour = "black"),
#axis.text.x = element_text(angle = 15,hjust = 1),
strip.text.x = element_text(size = 15,colour = "black"))+
geom_text(data = result.stepre.MF,aes(x=prd2,y=adjr2,label=adjsig),nudge_y = 0.05,size=7)+
geom_text(aes(label=round(adjr2,2)),vjust=2,color="black",size=3)+
facet_grid(~enzyme,scales = "free")
p3

p6<-ggplot(data=result.stepre.RF,aes(x=reorder(prd2, -adjr2),y=adjr2))+geom_col(fill="#fdae61",color=NA,width = 0.5)+
labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
axis.title = element_text(size = 15,colour = "black"),
#axis.text.x = element_text(angle = 15,hjust = 1),
strip.text.x = element_text(size = 15,colour = "black"))+
geom_text(data = result.stepre.RF,aes(x=prd2,y=adjr2,label=adjsig),nudge_y = 0.05,size=7)+
geom_text(aes(label=round(adjr2,2)),vjust=2,color="black",size=3)+
facet_grid(~enzyme,scales = "free")
p6

p9<-ggplot(data=result.stepre.DF,aes(x=reorder(prd2, -adjr2),y=adjr2))+geom_col(fill="#fdae61",color=NA,width = 0.5)+
labs(title = NULL, x = "Predictors", y = 'Adj. R-Square', fill = NULL)+
theme_bw()+theme(panel.grid = element_blank(),axis.text = element_text(size = 15,colour = "black"),
axis.title = element_text(size = 15,colour = "black"),
#axis.text.x = element_text(angle = 15,hjust =1),
strip.text.x = element_text(size = 15,colour = "black"))+
geom_text(data = result.stepre.DF,aes(x=prd2,y=adjr2,label=adjsig),nudge_y = 0.05,size=7)+
geom_text(aes(label=round(adjr2,2)),vjust=2,color="black",size=3)+
facet_grid(~enzyme,scales = "free")
p9

library(ggpubr)
ggarrange(p3,p6,p9,ncol = 3,nrow = 1,widths=c(1,1),heights=c(1,1),labels = c("(b)","(c)","(d)"))

#ggsave("plot2/Fig.5_stepwiseForest2.pdf",plot = last_plot(),units = "in",height = 3,width = 12,dpi = 600)
#Fig. S1
#SARSmbc#
g9<-ggboxplot(mydata.log,x='Treatment',y='SARSmbc',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific ARS activity\nLog(nmol ug MBC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g9

#SBG
g14<-ggboxplot(mydata,x='Treatment',y='SBG',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific BG activity\n(nmol ug MBC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g14

#SXYL
g15<-ggboxplot(mydata,x='Treatment',y='SXYL',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific XYL activity\n(nmol ug MBC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g15

#SNAG
g16<-ggboxplot(mydata,x='Treatment',y='SNAG',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific NAG activity\n(nmol ug MBC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g16

#SLAP
g17<-ggboxplot(mydata,x='Treatment',y='SLAP',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific LAP activity\n(nmol ug MBC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g17

#SAP
g18<-ggboxplot(mydata.log,x='Treatment',y='SACP',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific AP activity\nLog(nmol ug MBC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g18

ggarrange(g9,g14,g15,g16,g17,g18,ncol = 2,nrow = 3,widths=c(1,1),heights=c(1,1),
common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)","(e)","(f)"))

#ggsave("plot2/Fig.3_SpecificEnzymeActivity.pdf",plot = last_plot(),units = "in",height = 15,width = 14,dpi = 600)
#Fig. S2
#SARSt
g13<-ggboxplot(mydata,x='Treatment',y='SARSt',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific ARS activity\n(nmol ug TC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g13

#SBGt
mydata$SBGt<-mydata$BG/mydata$Total_C
g34<-ggboxplot(mydata,x='Treatment',y='SBGt',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific BG activity\n(nmol ug TC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g34

#SXYL
mydata$SXYLt<-mydata$XYL/mydata$Total_C
g35<-ggboxplot(mydata,x='Treatment',y='SXYLt',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific XYL activity\n(nmol ug TC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g35

#SNAG
mydata$SNAGt<-mydata$NAG/mydata$Total_C
mydata.log$SNAGt<-log(mydata$SNAGt+1)
g36<-ggboxplot(mydata.log,x='Treatment',y='SNAGt',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific NAG activity\nLog(nmol ug TC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g36

#SLAP
mydata$SLAPt<-mydata$LAP/mydata$Total_C
mydata.log$SLAPt<-log(mydata$SLAPt+1)
g37<-ggboxplot(mydata.log,x='Treatment',y='SLAPt',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific LAP activity\nLog(nmol ug TC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g37

#SAP
mydata$SACPt<-mydata$ACP/mydata$Total_C
mydata.log$SACPt<-log(mydata$SACPt+1)
g38<-ggboxplot(mydata.log,x='Treatment',y='SACPt',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Specific AP activity\nLog(nmol ug TC-1 h-1)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g38

ggarrange(g13,g34,g35,g36,g37,g38,ncol = 2,nrow = 3,widths=c(1,1),heights=c(1,1),
common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)","(e)","(f)"))

#ggsave("plot2/Fig.S2_SpecificEnzymeActivity2.pdf",plot = last_plot(),units = "in",height = 15,width = 14,dpi = 600)
#Fig. S3
#logS_C#
g10<-ggboxplot(mydata,x='Treatment',y='logS_C',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Log(ARS)/Log(BG+XYL)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g10

#logS_N#
g11<-ggboxplot(mydata,x='Treatment',y='logS_N',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Log(ARS)/Log(NAG+LAP)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g11

#logS_P#
g12<-ggboxplot(mydata,x='Treatment',y='logS_P',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Log(ARS)/Log(AP)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g12

#logC_N
mydata$logC_N<-log(mydata$BG+mydata$XYL)/log(mydata$NAG+mydata$LAP)
g31<-ggboxplot(mydata,x='Treatment',y='logC_N',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Log(BG+XYL)/Log(NAG+LAP)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g31

#logC_P
mydata$logC_P<-log(mydata$BG+mydata$XYL)/log(mydata$ACP)
g32<-ggboxplot(mydata,x='Treatment',y='logC_P',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Log(BG+XYL)/Log(AP)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g32

#logN_P
mydata$logN_P<-log(mydata$NAG+mydata$LAP)/log(mydata$ACP)
g33<-ggboxplot(mydata,x='Treatment',y='logN_P',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Log(NAG+LAP)/Log(AP)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g33

ggarrange(g10,g11,g12,g31,g32,g33,ncol = 2,nrow = 3,widths=c(1,1),heights=c(1,1),
common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)","(e)","(f)"))

#ggsave("plot2/Fig.4_enzymeActivity5.pdf",plot = last_plot(),units = "in",height = 15,width = 14,dpi = 600)
#Fig. S4
#TC
g19<-ggboxplot(mydata,x='Treatment',y='Total_C',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Total C (mg/g)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g19

#TN
g20<-ggboxplot(mydata,x='Treatment',y='Total_N',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Total N (mg/g)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g20

#TC_TN
g21<-ggboxplot(mydata,x='Treatment',y='TC_TN',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Total C/Total N")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g21

#WSOC
g22<-ggboxplot(mydata,x='Treatment',y='WSOC',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("WSOC (mg/kg)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g22

#NH4
g23<-ggboxplot(mydata,x='Treatment',y='NH4',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("NH4+-N (mg/kg)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g23

#NO3
g24<-ggboxplot(mydata,x='Treatment',y='NO3',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("NO3--N (mg/kg)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g24

ggarrange(g19,g20,g21,g22,g23,g24,ncol = 2,nrow = 3,widths=c(1,1),heights=c(1,1),
common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)","(e)","(f)"))

#ggsave("plot2/Fig.S3_soilproperties.pdf",plot = last_plot(),units = "in",height = 15,width = 14,dpi = 600)
#Fig. S5
#availP
g25<-ggboxplot(mydata.log,x='Treatment',y='availP',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Extractable P Log(mg/kg)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g25

#pH
g26<-ggboxplot(mydata,x='Treatment',y='pH',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("pH")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g26

#EC
g27<-ggboxplot(mydata.log,x='Treatment',y='EC',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("EC Log(us/cm)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g27

#MBC
#mydata$MBC<-mydata$MBC*1000
g28<-ggboxplot(mydata,x='Treatment',y='MBC',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("MBC (mg/g)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g28

#MBN
#mydata$MBN<-mydata$MBN/1000
g29<-ggboxplot(mydata,x='Treatment',y='MBN',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("MBN (mg/g)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g29

#MBC_MBN
g30<-ggboxplot(mydata.log,x='Treatment',y='MBC_MBN',fill = 'Treatment',facet.by = 'Forest',
add = "jitter",bxp.errorbar = T,bxp.errorbar.width = 0.3)+
scale_fill_manual(values = c("#31a354","#d8b365","#e66101","#b2abd2"))+
stat_compare_means(method = "t.test",comparisons = my_comparison)+
ylab("Log(MBC/MBN)")+
theme(axis.text = element_text(size = size,color = "black"),axis.title = element_text(size = size,color = "black"),
legend.text = element_text(size = size),legend.title = element_text(size = size),
strip.text.x = element_text(size = size))
g30

ggarrange(g25,g26,g27,g28,g29,g30,ncol = 2,nrow = 3,widths=c(1,1),heights=c(1,1),
common.legend = TRUE,labels = c("(a)","(b)","(c)","(d)","(e)","(f)"))

#ggsave("plot2/Fig.S4_soilproperties.pdf",plot = last_plot(),units = "in",height = 15,width = 14,dpi = 600)
#Fig. S6
ggplot(cor.out.all, aes(Y,X)) +
geom_tile(aes(fill = r), size=1)+
geom_text(data= cor.out.all[cor.out.all$r!=0,], aes(label = r), size = 5, font="bold")+
scale_fill_gradient2(guide = "legend", high='#2c7bb6',mid="white", low='#d7191c')+
theme_bw()+facet_grid(~Forest)+guides(fill = guide_legend(title = 'rho'))+
theme(axis.title= element_blank(),
axis.text.x=element_text(colour="black", size=15, face="bold",angle = 30, hjust=1),
axis.text.y=element_text(colour="black", size=15, face="bold"))

#ggsave("plot2/Fig.S1_spearmanCorr.pdf",plot = last_plot(),units = "in",height = 6,width = 12,dpi = 600)
R Markdown